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Abstract 

Cascading failures in interdependent networks have been investigated using 
percolation theory in recent years. Here, we study the dynamics of the cascad- 
ing failures, the average and fluctuations of the number of cascading as a function 
of system size N near criticality. The system we analyzed is a pair of fully in- 
terdependent Erdos-Renyi (ER) networks. We show that when p is close to p c , 
the whole dynamical process of cascading failures can be divided into three time 
stages. The giant component sizes in the second time stage, presented by a plateau 
in the size of giant component, have large standard deviations, which cannot be 
well predicted by the mean-field theory. We also investigate the standard deviation 
of the total time std(r) using simulations. When p = p c , our numerical simulations 
indicate that std(r) ~ A' 1 ' 3 , which increases faster than the mean, < r >~ N ,/4 , 
predicted by the mean-field theory. We also find the scaling behavior as a function 
of N and p of < r > and std(r) for p < p c . 

1 Introduction 

It has been shown that many complex systems in the real world are interdependent 
and interact with each other 1 1 2,3 4], The cascading process of failures in two 
fully interdependent networks has been investigated by Buldyrev et al. using both 
theory and simulations [5|. It has been shown that when a fraction 1 - p of nodes 
are randomly removed in one network, assuming that only the giant component 
of each network still functions, it yields cascading failures due to interdependency 
and percolation processes. The size of the final mutual giant component collapses 
abruptly as in a first order phase transition when reducing p [5] [6). Leicht and 
D'Souza also provided a theoretical framework to study the percolation in a system 
of interacting networks (7). 
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Several further studies on percolation of interdependent networks have been 
performed. For example, Parshani et al. have investigated the effects of the cou- 
pling strength q between the networks (the fraction of interdependent nodes in both 
networks) on the phase transition properties (8). It is found that there is a critical 
dependency q c above which the transition is first order and below it is continuous 
as in second order phase transitions 1 8 1. Shao et al. have extended the theoretical 
analysis to multiple support-dependent links (5). Hu et al. have studied the case 
where both interconnectivity links and interdependency links exist 1 10 1. Gao et al. 
generalized the study of a pair of interdependent networks to a network of coupled 
networks (NON) t6l llll . Huang et al. and Gao et al. obtained analytical results 
for fully and partially interdependent networks under target attacks, respectively 
1 1211 1 31 . Very recently, percolation on interdependent lattices has been also inves- 
tigated with the finding that the vulnerability of spatially embedded interdependent 
networks is significantly higher compared to random networks 1 1411 151 . 

In these studies, researchers mainly focus on theoretical results based on the 
mean-field (MF) approach, which assumes that the number of nodes N — > oo, and 
compare them with simulations on systems of size N. Buldyrev et al. analysed the 
effects of N on the average total time (the total number of cascades) |5|. In this 
paper, we mainly discuss the stages in the cascading process, and study the varia- 
tions of the total time to collapse for networks of size A' for fully interdependent 
Erdos-Renyi (ER) networks (q = 1). 

2 Three Stages in the Dynamical Process 

The model of fully interdependent ER networks is defined as follow. Assume that 
A and B are two ER networks of the same size A'. The average degrees of A and B 
are k A and k B . Each A-node a, depends on exactly one randomly-chosen S-node bj, 
and bj also only depends on a,. The initial attack is randomly removing a fraction 
1 - p of A-nodes. As nodes and edges are removed, each network breaks up into 
connected components (clusters). Assume that when the network is fragmented, 
the nodes belonging to the largest component (giant component) connecting a fi- 
nite fraction of the network are still functional. Since each network is connected 
differently, the nodes that become nonfunctional on each step are different for both 
networks. This will cause a cascading process of failures in the system. 

First, we investigate the variations of the giant component size tf/, of network 
A at time step t during the dynamical cascading failure process. Fig. |l(a)| shows 
the giant component of network A in both theory and simulation when p is close 
to but below p c . 

When considering fifteen realizations in simulations in Fig. |l(a)| we can see 
that they fall into two types according to whether the system totally collapses ul- 
timately. Only in six realizations, the system fully collapses, and the total time r 
for each of them to collapse is not the same. The dynamical process for these 6 
realizations can be divided into three stages as follow: in Stage 1, the giant compo- 
nent decreases relatively fast after the initial attack; in Stage 2 (a plateau), the giant 
component size decreases very slowly; in Stage 3, the system goes toward a total 
collapse relatively fast. For the other nine realizations where the system does not 
totally collapse, we can also observe a similar Stage 1, and the dynamical process 
ends in Stage 2 with a non-zero mutual giant component. 

Fig. |l(b)| illustrates how the standard deviation of the giant component size iff, 
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Figure 1 : a. Dynamical process of the giant component size ipt for fully coupled ER 
networks in both theory and simulation (15 realizations), k — 5, p — 0.4908 for theory, 
p = 0.49108 for simulation, N = 250,000. b. Mean and standard deviation of the 
giant component size in the first several time steps for coupled ER networks, k = 5, 
p = 0.49108, N = 500,000, M = 10,000. 

varies with t in the first time steps for M = 10, 000 realizations. 

While in Stage 1 we find very small variations, in Stage 2, we find large varia- 
tions of the giant component sizes between different realizations. This means that 
the MF theory cannot well predict how long in time (how many iterations of the 
cascading failures) appear in Stage 2. 

3 Variations of r 

As discussed in the previous section, the fluctuations of Stage 2 cannot be predicted 
by the MF theory. This is the feature of the fluctuations of the total time t, at 
which different realizations of the system completely collapse. Here, we mainly 
investigate the role of p and N on the standard deviation of r. Buldyrev et al. 
analyzed the scaling behavior of the average total time < r > |5|. For systems 
of size N with p < p c , it was shown that the mean of r diverges at p c as < t >~ 
1/ sjPc - P\ for p = p c , < t >~ N l/4 ; for p > p c , < r >~ InN/ y/p - p c . 

We wish to consider here the fluctuations in r. In our simulation, we consider 
the case where p < p c , and only those realizations which totally collapses. We wish 
to better understand how N and p influence the mean and the standard deviation of 
the total time t of the totally collapsed systems. 

Fig. |2(a) | and |2(b) | show the effects of N and p (near p c ) on the mean and the 
standard deviation of r in simulation respectively. It can be seen from Fig. |2(a)| 
that < t > increases with N as N ~ N 1 ' 4 when p = p c ~ 0.491082 |5|. However, 
when p < p c , < t > becomes a constant for larger values of N. According to these 
behaviors, we assume the following scaling function, 

< t >~ N 1 ' 4 ■ f(u). (1) 

Here, u = (p c - p) ■ N [/a , and /(«) is a function of u that satisfies: f(u) is a constant 
for u « 1, and f(u) ~ iT" 14 for u » 1. 

Fig. |2(c)| is a scaled version of Fig. |2(a)| by plotting < r > /TV 1 ' 4 versus 
u = (p c - p) ■ N l!a , which aims to test the assumption of the scaling form of Eq. [TJ 
Here, two more values of p close to p c are included: p = 0.4908 and p = 0.491. 
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We find that the best choice of a for obtaining a good scaling is a = 2. Thus, Fig. 
|2(c) | supports the scaling form of Eq. (1). In this way, we can see that the slope of 
each curve changes from to about -1/2 at u = (p c - p) ■ N [/1 ~ 1. Therefore, the 
scaling behavior of < r > is for N > N* ~ (p c - p)~ 2 , 

< r >~ N 114 ■ u 112 = (p c - p)- 1 ' 2 , (2) 
independent of N (Fig. |2(a)fr , and for N < N*, 

<t>~ N l/ \ (3) 

independent of p (Fig. |2(a)[ >. Note that Eqs. [2] and [3] are consistent with those of 
Buldyrev et al. (5J. Here, we suggest the full scaling behavior, Eq. Q~]which yields 
the crossover N", 

AT ~ ( Pc - p)-" = (p c - p y\ (4) 

between the critical behavior for N < N" and non-critical for N > N" . As seen for 
P ~ * Pa N* — > oo and for all N one observe a critical behavior. 

Fig. |2(b)| illustrates the effects of N and p on the standard deviation std(r). 
First, our numerical simulations suggest that for p = p c , the slope of std(r) with N 
is about 1/3, which means std(r) ~ A' 1 ' 3 , i.e., it increases faster than the mean, at 
least in this range of N. Furthermore, when p < p c , std(r) converges for larger N, 
and the slope is about -2/3. Similarly to < t >, we assume the scaling form 

std(r) ~ A' 1 ' 3 • g(u), (5) 

where u = (p c — p) • N [/a , and g(u) satisfies: g{u) is a constant for u « 1, and 
g(u) ~ u~ a for u » 1. 

Fig. |2(d)| shows the scaling behavior of std(r) assumed in Eq. |5]is supported 
by simulations with the best choice of a again about 2. The slope of the right tail 
in Fig. |2(d)| is about -2. Therefore, for N > N* ~ (p - p c Y 2 , we have 

std(r) ~ N [p ■ u 1 = (p c - pT 2 ■ AT 2 ' 3 , (6) 

and for N < N", 

std(r) ~ N U3 , (7) 

both consistent with Fig. |2(b)| 

Our results indicate that although the value of p c in simulation converges to 
the theoretical value, the total time t caused by complete collapses at p = p c 
does not converge, and varies from one realization to another. In order to better 
understand this behavior, we performed the following test on the total time r in 
theory. First, we obtain the distribution of p c (denoted by p s c to avoid confusion 
with p c in theory) in simulations for each value of A'. Fig. [3]shows the mean and 
the standard deviation of pf, in 3, 000 realizations for different values of Af. Notice 
that here, the number of initially randomly removed nodes is fixed at p = p c . We 
find that < pf, > is almost a constant, and std(pf ) ~ A' -049 . This result reflects 
the randomness that results from the random initial attack as well as the random 
network structure. By linear fitting on both < p\ > and std(p^ ) respectively, we 
can predict the distribution of p s c assuming a normal distribution for any system 
size N. Next, we artificially use this distribution of p s c as the distribution of p in 
the initial attack for the theoretical dynamical process (Eq. 7 in SI of |5 1). In this 



4 




Figure 2: a. Effects of the system size N on the average total time r for different values 
of p. k = 5, M — 3,000. b. Effects of the system size N on the standard deviation 
of the total time r for different values of p. k = 5, M — 3,000. c. Scaling behavior 
of the average r. Two more values of p are included: p = 0.4908 and p = 0.491. d. 
Scaling behavior of the standard deviation of r. Two more values of p are included: 
p = 0.4908 and p = 0.491. 



Figure 3: (Color online) Mean and standard deviation of p c in simulation, k = 5, 
M = 3, 000. The number of initially attacked nodes is fixed. 
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Figure 4: (Color online) Mean and standard deviation of r in theory, k — 5, 
M = 10,000. For all values of N, use linear fitting on the mean and the standard 
deviation of p s c obtained in simulation, respectively, and p is chosen according to a 
normal distribution using these extended means and standard deviations. 



way, we introduce randomness into the theory. Fig. [4] shows that in this modified 
theoretical process, both < r > and std(r) increase as close to N >14 . Compared 
with our simulation results (Fig. |2(b)[ (, this behavior indicates that apart from the 
random initial attack, other sources that cause more randomness on the total time 
should exist, like the network topologies and the critical behavior in Stage 2. 
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